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ABSTRACT 

On megaparsec scales the Universe is permeated by an intricate filigree of clusters, 
filaments, sheets and voids, the Cosmic Web. For the understanding of its dynamical 
and hierarchical history it is crucial to identify objectively its complex morphological 
components. One of the most characteristic aspects is that of the dominant underdense 
Voids, the product of a hierarchical process driven by the collapse of minor voids in 
addition to the merging of large ones. 

In this study we present an objective void finder technique which involves a min- 
imum of assumptions about the scale, structure and shape of voids. Our void finding 
method, the Watershed Void Finder (WVF), is based upon the Watershed Transform, 
a well-known technique for the segmentation of images. Importantly, the technique 
has the potential to trace the existing manifestations of a void hierarchy. The basic 
watershed transform is augmented by a variety of correction procedures to remove 
spurious structure resulting from sampling noise. 

This study contains a detailed description of the WVF. We demonstrate how it 
is able to trace and identify, relatively parameter free, voids and their surrounding 
(filamentary and planar) boundaries. We test the technique on a set of Kinematic 
Voronoi models, heuristic spatial models for a cellular distribution of matter. Com- 
parison of the WVF segmentations of low noise and high noise Voronoi models with 
the quantitatively known spatial characteristics of the intrinsic Voronoi tessellation 
shows that the size and shape of the voids are succesfully retrieved. WVF manages to 
even reproduce the full void size distribution function. 

Key words: Cosmology: theory - large-scale structure of Universe - Methods: data 
analysis - numerical 
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1 INTRODUCTION 

Voids form a prominent aspect of the distribution of galax- 
ies and matter on megaparsec scales. They are enormous 
regions with sizes in the range of 20 — 50ft" 1 Mpc that 
are practically devoid of any galaxy and usually roundish in 
shape. Forming a n essential ingredient of the Cosmic Web 
|Bond et fflZj|l99q ). they are surrounded by elongated fila- 
ments, sheetlike walls and dense compact clusters. Together 
they define the salient weblike pattern of galaxies and mat- 
ter which pervades the observable Universe. 

Voids have been known as a featur e of galaxy surveys 
since the first surveys were compiled ( Chincar ni fc Roodl 
ll975l : lGregorv fc Thompsonll 19781: lEinasto et ql.ll 19801) . Fol- 
lowing the discovery by ( Kirshner et qZ.Tl98ll . 19871 ) of the 
most dramatic specimen, the Bootes void, a hint of their 
central position within a w eblike arrangement came with 
the f irst CfA redshift slice l|de Lapparent Geller fc Huchral 
1986). This view has recently been expanded dramatically 
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as maps of the spatial distrib ution of hundreds o f thousands 
of galaxies in the 2dF GRS l|Colless et aZj|2003h and SDSS 
ijAbazaiian et al1l2003l ) have become available. 

Voids are a manifestation of the cosmic structure 
formation process as it reaches a non-linear stage of 
evolution. Structure forms by gravitational instability from 
a primordial Gaussian field of small amplitude density 
pertu rbati o ns, w h ere voids emerge out of the dep r ession s 
(e.g. [fcki Il984l : Ivan de Weveaert fc van Kampenl Il993l ). 
They mark the transition scale at which perturbations have 
decoupled from the Hubble flow and organized themselves 
into recognizable structur al features. Early theor e tical 
mode l s of void formation llHoffman fc Shaham Il982l : llckel 
1 19841 : lBertschinge"rl 1 19851 : iBlumenthal et ali 1 199*3 ) were 
followed and generalized by the fi rst numerical s imula- 

Gellerl 



tions of void centered universes 



van dc Wcygacrt & van Kampen 



Regos fc Gellerl 119911 
199*"l : iDubinski et all 



1993t iMartel fc Wassermannlll990l ). 

In recent years the huge increase in computa- 
tional resources has enabled N-body simulations to re- 
solve in detail the intricate substructure of voids within 
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the context of h ierarchical cosmologica l structure for- 

mation scenarios dMathis fc White! 2002 ; Gottlober et all 
■ i 1*-^ — ^ — 



20031: [Hoeft et al 



"[20071: lArbabi-Bidgoli fc Mullerl 120021 
Goldberg fc Vogelevll20ol;IColberg et a/j|2005l;|Padirla et all 



2005). They confirm the theoretical expectation of voids 



having a rich substructure as a result of their hierarchi- 
cal buildup. Theoretically this evolution has been succes- 
fully embedded in the extended Press-Schechter descrip - 
tion ilPress fc Schechterll974l;lBond et ffll.ll99ll ;I Shethll998n . 
ISheth fc van de WeygaertJ ( 2004 ) showed how this can be 



des cribed by a two-barrier excursion set formalism (also 
see iFurlanetto fc Pir an 2006). The two barriers refer to the 
two processes dictating the evolution of voids: their merg- 
ing into ever larger voids as well as the collapse and disap- 
pearance of small ones embe dded in overdense regions (see 
Ivan de Weygaert et ai\\2004 ). 

Besides representing a key constituent of the cos- 
mic matter distribution voids are interesting and impor- 
tant for a variety of reasons. First, they are a promi- 
nent feature of the megaparsec Universe. A proper and 
full understanding of the formation and dynamics of the 
Cosmic Web is not possible w ithout understanding the 
struc ture and evolution of voids l|Sheth fc van de Weygaert] 
120041 '). Secondly, they are a probe of cosmological param- 
eters. The outflow from the voids depends on the mat- 
ter density parameter f2 m , the Hubble parameter H(t) 



and possibly on the cosmoloj 


4cal constant A (see e.g. 


Bernardeau & van de Weygaert 


19961: Dekel & Reesl Il994l; 


Martcl & Wasscrmann 1990: Flichc & Triay 2006). These 



parameters also dicta t e their redshift sp ace distortions 
jRvden fc Melott]|l996l ; ISchmidt e£~oZll200 j ) while their in- 
trinsic structure and shape is sensitive to var ious aspects 
of th e power spectrum of density fluctuations l|Lee fc Park! 
2006). A third point of interest concerns the galaxies 
in voids. Voids provide a unique and still largely pris- 
tine environment for s t udying the evolution of galax- 
ies llHoffman et all Il992l : iLittle fc Weinberg! Il994l : iPeebles! 
l200lT T The recent interest in environmental influences 



in this direction (Szomoru et al. 




19981; Grogin & Gellerl 


19991; iMathis & White 20021; Friedmann & Piranl 200ll; 


Benson et all 1996; Gottlober et al. 20031; Hoeft et all 


2007; Furlanetto & Piran 2006; 


Hovle & Vogclcy 2002; 


Roias et aZ.ll2005l: iPatiri et all 2006; 


Ceccarelli et aZ.I:2006). 



Despite the considerable interest in voids a fairly ba- 
sic yet highly significant issue remains: identifying voids 
and tracing their outline within the complex spatial ge- 
ometry of the Cosmic Web. There is not an unequivo- 
cal definition of what a void is and as a result there 
is considerable disagree ment on the precise outline of 
such a region (see e.g. IShandarin et all [2006). Because 
of the vague and diverse definitions, and the diverse in- 
terests in voids , there is a plethora of void identifica- 
tion pr ocedures llKauffmann fc Fairall 1991 ; El- Ad fc Piranl 



19971; I Aikio fc Mahonenl 119981; iHovle fc Vogelevl 1 2002 



Arba bi-Bidgoli & MulleT l200a IPlionis fc Basilakosl 
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Patiri et aZ.ll2006l ; ICoiberg et aZ.l2005l ; IShandarin et aZj200e 



2002 



Hahn et aZj|2007l : lNevrinckll2007h . 

The "sp h ere-ba sed" voidfinder algorithm of 
lEl-Ad fc Piranl (|l997h has been at the basis of most 
voidfinding methods. However, this succesful approach 
will not be able to analyze complex spatial configurations 



in which voids may have arbitrary shapes and contain a 
range and variety of substructures. A somewhat related 
and tessellation based voidfindi ng technique t hat still is 
under development is ZOBOV (iNevrinckl |2007h . It is the 
voidfinder equivalent to the VOBO Z halofinder method 
jNevrinck, Gnedin fc Hamilton! I2005T ). 

Here we introduce and test a new and objective 
voidfinding formalism that has been specifically designed 
to dissect the multiscale character of the void network and 
the weblike features marking its boundaries. Our Watershed 
Void Finder (WVF) is based on the watershed algo rithm 
dBeucher fc Lantueioull Il979l ; iBeucher fc Meverl Il993l ). It 
stems from the field of mathematical morphology and image 
analysis. 

The WVF is defined with respect to the 
DTFE density field of a dis crete point distribution 
( Sch aap fc van de Wevgaertll200Ch . This assures an optimal 
sensitivity to the morphology of spatial structures and 
yields an unbiase d probe of substructure in the mass dis- 
tribu tion (see e.g lOkabel |2000| ; ISchaap fc van de Weygaert] 
2000). Because the WVF void finder does not impose a 
priori constraints on the size, morphology and shape of 
a voids it provides a basis for analyzing the intricacies of 
an evolving void hierarchy. Indeed, this has been a major 
incentive towards its development. 

This study is the first in a series. Here we will define 
and describe the Watershed Void Finder and investigate its 
performance with respect to a test model of spatial we- 
blike distributions, Voronoi kinematic models. Having as- 
sured the success of WVF to trace and measure the spatial 
characteristics of these models the follow-up study will ad- 
dress the application of WVF on a number of GIF N-bod y 
simulations of structure formation l|Kauffmann et allll999l ). 
Amongst others, WVF will be directed towards character- 
izing the hierarchical structure of the me gaparsec void pop- 
ulation |Sheth fc van de Wevgaert]|2004r i. For a comparison 
of the WVF with other void finder meth o ds we refer to the 
extensive study of IColberg. Pearce et al] (|2007h - 

In the following sections we will first describe how the 
fundamental concepts of mathematical morphology have 
been translated into a tool for the analysis of cosmologi- 
cal density fields inferred from a discrete N-body simulation 
or galaxy redshift survey point distribution (sect. 2 & 3). To 
test our method we have applied it to a set of heuristic and 
flexible models of a cellular spatial distribution of points, 
Voronoi clustering models. These are described in section 4. 
In section 5 we present the quantitative analysis of our test 
results and a comparison with the known intrinsic proper- 
ties of the test models. In section 6 we evaluate our findings 
and discuss the prospects for the analysis of cosmological 
N-body simulations. 



2 THE WATERSHED VOID FINDER 

The new void finding algorithm whi ch we introduce here is 
based on t he watershed transform o f IBeucher fc Lantueioul 
(1979) and IBeucher fc Meverl (|l993h . A more extensive and 
technical description of the basic concepts of mathematical 
morphology and the basic watershed algorithm in t erms of 
homotopy transformations on lattices (Krcsch 1998) is pro- 
vided in appendix |A~1 and |B1 
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Figure 1. Three frames illustrating the principle of the watershed transform. The lefthand frame shows the surface to be segmented. 
Starting from the local minima the surrounding basins of the surface start to flood as the water level continues to rise (dotted plane 
initially below the surface). Where two basins meet up near a ridge of the density surface, a "dam" is erected (central frame). Ultimately, 
the entire surface is flooded, leaving a network of dams defines a segmented volume and delineates the corresponding cosmic web 
(righthand frame). 



2.1 the Watershed Transform (WST) 

The watershed transform is used for segmenting images 
into distinct regions and objects. The Watershed Trans- 
form (WST) is a concept defined within the context of 
mathematical morpho l ogy, and was first introduced by 
iBeucher fc Lantueioull (| 19791 ). The basic idea behind the 
WST finds its origin in geophysics. The WST delineates 
the boundaries of the separate domains, i.e. the basins, into 
which yields of, for example, rainfall will collect. 

The word watershed refers to the analogy of a landscape 
being flooded by a rising level of water. Suppose we have a 
surface in the shape of a landscape (first image of Fig. [JJ. 
The surface is pierced at the location of each of the min- 
ima. As the water-level rises a growing fraction of the land- 
scape will be flooded by the water in the expanding basins. 
Ultimately basins will meet at the ridges corresponding to 
saddle-points in the density field. This intermediate step is 
plotted in the second image of Fig. [JJ The ridges define the 
boundaries of the basins, enforced by means of a sufficiently 
high dam. The final result (see last in Fig. [TJ of the com- 
pletely immersed landscape is a division of the landscape 
into individual cells, separated by the ridge dams. In the 
remainder of this study we will use the word "segment" to 
describe the watershed's cells. 



2.2 Watershed segments: qualities 

The watershed algorithm holds several advantages with re- 
spect to other voidfinders: 

• Within an ideal smooth density field (i.e. without noise) 
it will identify voids in a parameter free way. No predefined 
values have to be introduced. In less ideal, and realistic, cir- 
cumstances a few parameters have to be set for filtering out 
discreteness noise. Their values are guided by the properties 
of the data. 

• The watershed works directly on the topology of the 
field and does not reply on a predefined geometry /shape. 
By implication the identified voids may have any shape. 



• The watershed naturally places the divide lines on the 
crests of a field. The void boundary will be detected even 
when its boundary is distorted. 

• The transform naturally produces closed contours. As 
long as minima are well chosen the watershed transform will 
not be sensitive to local protrusions between two adjacent 
voids. 

Obviously we can only extract structural information to the 
extent that the point distribution reflects the underlying 
structure. Undersampling and shotnoise always conspire to 
obfiscate the results, but we believe the present methodology 
provides an excellent way of handling this. 



2.3 Voids and watersheds 

The Watershed Void Finder (WVF) is an implementation 
of the watershed transform within a cosmological context. 
The watershed method is perfectly suited to study the holes 
and boundaries in the distribution of galaxies, and holds the 
specific promise of being able to recognize the void hierarchy 
that has been the incentive for our study. 

The analogy of the WST with the cosmological context 
is straightforward: voids are to be identified with the basins, 
while the filaments and walls of the cosmic web are the ridges 
separating the voids from each other. 



2.4 The Watershed Void Finder: Outline 

An outline of the steps of the watershed procedure within 
its cosmological context is as follows: 

• DTFE: Given a point distribution (N-body, redshift 
survey), the Delaunay Tessellatio n Field Estimator (DTFE, 
ISchaap fc van de Wevgaertj |2000| ) is used to define a con- 
tinuous density field throughout the sample volume. This 
guarantees a density field which retains the morphological 
character of the underlying point distribution, i.e. the 
hierarchical nature, the web-like morphology dominated by 
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filaments and walls, and the presence voids is warranted. 

• Grid Sampling: For practical processing purposes 
the DTFE field is sampled on a grid. The optimal grid size 
has to assure the resolution of all morphological structures 
while minimizing the number of needed gridcells. This 
criterion suggests a grid with gridcells whose size is in the 
order of the interparticle separation. 

• Rank-Ordered Filtering: The DTFE density field 
is adaptively smoothed by means of Natural Neighbour 
Maxmin and Median filtering. This involves the compu- 
tation of the median, minimum or maximum of densities 
within the contiguous Voronoi cell, the region defined by a 
point and its natural neighbours . 

• Contour Levels: The image is transformed into a 
discrete set of density levels. The levels are defined by a 
uniform partitioning of the cumulative density distribution. 

• Pixel Noise: With an opening and closing (operation 
to be defined in appendix. |X} of 2 pixel radius we further 
reduce pixel by pixel fluctuations. 

• Field Minima: The minima in the smoothed density 
field are identified as the pixels (grid cells) which are 
exclusively surrounded by neighbouring grid-cells with a 
higher density value. 

• Flooding: The flooding procedure starts at the lo- 
cation of the minima. At successively increasing flood 
levels the surrounding region with a density lower than 
the corresponding density threshold is added to the basin 
of a particular minimum. The flooding is illustrated in Fig.[T] 

• Segmentation: Once a pixel is reached by two distinct 
basins it is identified as belonging to their segmentation 
boundary. By continuing this procedure up to the maximum 
density level the whole region has been segmented into 
distinct void patches. 

• Hierarchy Correction: A correction is necessary to 
deal with effects related to the intrinsic hierarchical nature 
of the void distribution. The correction involves the removal 
of segmentation boundaries whose density is lower than 
some density threshold. The natural threshold value would 
be the typical void underdensity A = —0.8 (see sect. l3.47Tj) . 
Alternatively, dependent on the application, one may chose 
to take a user-defined value. 



2.5 WVF by example: 

Voids in a ACDM simulation 

A direct impression of the watershed voidfinding method is 
most readily obtained via the illustration of a representative 
example. In Fig. [5] the watershed proced ure has been ap- 
plied to the cosmological GIF2 simulation (|Kauffmann et all 
1 19991 ). 

The N-body particle distribution (lefthand Fig. [2J is 
translated into a density field using the DTFE method. 
The application of the DTFE method is described in sec- 



tion 13.11 the details of the DTFE procedure are specified in 
Appendix [Dl 

The DTFE density field is sampled and interpolated on 
a 256 3 grid, the result of which is shown in the top right- 
hand frame of Fig. [2] The gray-scales are fixed by uniformly 
sampling the cumulative density distribution, ensuring that 
all grayscale values have the same amount of volume. 

The DTFE density field is smoothed by means of the 
adaptive Natural Neighbour Median filtering described in 
sect. 13.21 This procedure determines the filtered density 
values at the location of the particles. Subsequently, these 
are interpolated onto a grid. This field is translated into a 
grayscale image following the same procedure as that for the 
raw DTFE image (bottom lefthand panel). 

The minima in the smoothed density field are identi- 
fied and marked as the flooding centres for the watershed 
transform. The resulting WVF segmentation is shown in the 
bottom righthand frame of Fig. [2] 

The correspondence between the Cosmic Web, its voids 
and the watershed segmentation is striking. There is an al- 
most perfect one-to-one correspondence between the seg- 
mentation and the void regions in the underlying density 
field. The WVF method does not depend on any predefined 
shape. As a result, the recovered voids do follow their natural 
shape. A qualitative assessment of the whole simulation cube 
reveals that voids are very elongated and have a preferential 
orientation within the cosmic web, perhaps dictated b y the 
megaparsec tidal force field (see e.g. lLee fc Parkl ;2006V 

Clearly, the Watershed Void Finder is able to extract 
substructure at any level present in the density distribution. 
While this is an advantage with respect to tracing the pres- 
ence of substructure within voids it does turn into a disad- 
vantage when seeking to trace the outline of large scale voids 
or when dealing with noise in the dataset. While the noise- 
induced artificial segments are suppresed by means of the 
full machinery of Markers (sect. 13. Void Patch Merging 
(sect. |3~4)) and Natural Neighbour Rank filtering (sect. [3~2l) . 
it are the latter two which may deal with intrinsic void hi- 
erarchy. 

The follow-up study (|Platen. van de Weygaert fc Jonesl 

l2007fl will involve a detailed quantitative analyze of volume 
and shapes of the voids in the GIF2 mass distribution for a 
sequence of timesteps. 



3 METHOD: DETAILED DESCRIPTION 

In order to appreciate the various steps of the Watershed 
Void Finder outlined in the previous section we need to de- 
scribe a few of the essential steps in more detail. 

To process a point sample into a spatial density field 
we use DTFE. To detect voids of a particular scale it is 
necessary to remove statistically insignificant voids gener- 
ated by the shotnoise of the discrete point sample as well as 
physically significant subvoids. In order to retain only the 
statistically significicant voids we introduce and apply Nat- 
ural Neighbour Rank-Order filtering. Hierarchy Merging is 
used for the removal of subvoids which one would wish to 
exclude from a specific void study. 
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Figure 2. A visualization of several intermediate steps of the Watershed VoidFinding method. The top lefthand frame shows the 
particles of a slice in the LCDM GIF simulation. The corresponding DTFE density field is shown in the top righthand frame. The next, 
bottom lefthand, frame shows the resulting 5th order median-filtered image. Bottom righthand frame: the resulting WVF segmentation, 
computed on the basis of the median filtered image. The image shows the superposition of WVF ridges (black) on the original density 
field. 



3.1 The DTFE density field 

The input samples for our analysis are mostly samples of 
galaxy positions obtained by galaxy redshift surveys or the 
positions of a large number of particles produced by N-body 
simulations of cosmic structure formation. In order to de- 
fine a proper continuous field from a discrete distribution 
of points - computer particles or galaxies - we translate 
the spatial point sample into a continuous density field by 



means of the Delaunay Tessellatio n Field Estimator (DTFE, 
Sch aap fc van de Wevgaertll2000h . 



.1.1 DTFE 



The DTFE technique l|Schaap fc van de Wevgaertll2000l ) re- 
covers fully volume-covering and volume-weighted continu- 
ous fields from a discrete set o f sampled field values. The 
method has been developed by ISchaap fc van de Weygaerj 
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Figure 3. Natural Neighbours of a point. The black dot rep- 
resents the central point, the open circles its Natural Neigh- 
bours. The solid edges mark the Voronoi cell surrounding the cen- 
tral point, along with the connecting Voronoi edges. The dashed 
lines delineate the corresponding Delaunay triangles. The central 
Voronoi cell is surrounded by its related Delaunay triangles, defin- 
ing the Natural Neighbours. The image is an illustration of the 
dual relationship between Voronoi and Delaunay tessellations. 



(2000) and forms an elabo ration of the velocity interpola- 
tion s cheme introduced by iBernardeau fc van de Weygaert! 
(1996). It is based upon the use of the Voronoi and Delau- 
nay tessellations of a given spatial point distribution to form 
the basis of a natural, fully self-adaptive filter in which the 
Delaunay tessellations are used as multidimensional inter- 
polation intervals. A typical example of a DTFE processed 
field is the one shown in the top r ow of Fig. [2] the parti- 
cles of a GIF N-body simulation l|Kauffmann et all Il999l ) 
are translated into the continuous density field in the right- 
hand frame. 

The primary ingredient of the DTFE method is the De- 
launay tessellation of the particle distribution. The Delau- 
nay tessellation of a point set is the uniquely defined and 
volume-covering tessellation of mutually disjunct Delaunay 
tetrahedra (triangles in 2D). Each is defined by the set of 
four points whose circumscribing sphere does n ot contain 
any of the other points in the generating set l|Delaunavl 
1 19341 ). The Delaunay tessellation and the Voronoi tessel- 
lation of the point set are each others dual. The Voronoi 
tessellation is the division of space into mutually disjunct 
polyhedra, each polyhedron consisting of the part of space 
closer to the defining poin t than any of the other points 
|Voronoilll90sl ; IOkabeil2000l ) 

DTFE exploi ts three prop e rties of Voronoi and Delau- 
nay t essellations (|Schaadl2007l ; ISchaap fc van de WeygaertJ 
120071 ). The tessellations are very sensitive to the local point 
density. DTFE uses this to define a local estimate of the 
density on the basis of the inverse of the volume of the tes- 
sellation cells. Equally important is their sensitivity to the 
local geometry of the point distribution. This allows them 
to trace anisotropic features such as encountered in the cos- 



Figure 4. Examples of 2-D grid connectivities. By default the 
central square is white. Cells connected to the centre are repre- 
sented by gray squares. Lefthand frame: a 4-connectivity. Centre 
frame: a 8-connectivity. Righthand frame: a structure element 
representing a ball of 2 pixels. 



mic web. Finally, DTFE exploits the adaptive and minimum 
triangulation properties of Delaunay tessellations in using 
them as adaptive spatial interpolation intervals for irregular 
point distributions. In this way i t is the first order version 
of the Natural Neighbour met hod l|Braun fc S ambridgel ll99"5l ; 
ISukumar|[l998l ; IWatsorJll992l ). 

Within the cosmological context a major - and crucial 
- characteristic of a processed DTFE density field is that it 
is capable of delineating three fundamental characteristics 
of the spatial structure of the megaparsec cosmic matter 
distribution. It outlines the full hierarchy of substructures 
present in the sampling point distribution, relating to the 
standard view of structure in the Universe having arisen 
through the gradual hierarchical buildup of matter concen- 
trations. DTFE also reproduces any anisotropic patterns in 
the density distribution without diluting their intrinsic ge- 
ometrical properties. This is particularly important when 
analyzing the the prominent filamentary and planar fea- 
tures marking the Cosmic Web. A third important aspect 
of DTFE is that it outlines the presence and shape of void- 
like regions. Because of the interpolation definition of the 
DTFE field reconstruction voids are rendered as regions of 
slowly varying and moderately low density values. 

A more detailed outline of the DTFE reconstruction 
procedure can be found in appendix [Dl 



3.1.2 DTFE grid 

DTFE involves the estimate of a continuous field through- 
out the complete sample volume. To process the DTFE field 
through the WVF machinery we sample the field on a grid. 
It is important to choose a grid which is optimally suited for 
the void finding purpose of the WVF method. On the one 
hand, the grid values should represent all physically signif- 
icant structural features (voids) in the sample volume. On 
the other hand, the grid needs to be as coarse as possible 
in order to suppress the detection of spurious and insignifi- 
cant features. The latter is also beneficial from a viewpoint 
of computational efficiency. This is achieved by adopting a 
gridsize in the order of the mean interparticle distance. 

The DTFE grid sampling is accomplished through 
Monte Carlo sampling within each grid cell. Within each 
gridcell the DTFE density value is measured at 10 randomly 
distributed sample points. The grid value is taken to be their 
average. 
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Figure 5. The concept of hierarchical watershed. Not all divide lines produced by the watershed may be relevant. They are removed if 
they do not fulfil a particular criterium (e.g. if they have a contrast lower than some threshold). Only the significant watershed segments 
survive. The segmentation after 5 iterative density smoothings and removal of boundaries below a contrast of 0.8. 



3.2 Natural Neighbour Rank-Ordered Altering 

A major and novel ingredient of our WVF method intended 
to eliminate shot noise in the DTFE density field reconstruc- 
tions is that of a natural non-linear filtering extension: the 
Natural Neighbour Rank- Ordered filtering 

We invoke two kinds of non-linear adaptive smoothing 
techniques, Median Filtering and Max/Min Filtering, the 
latter originating in mathematical morphology (MM). Both 
filters are rank order filters, and both have well known be- 
haviour. They have a few important properties relevant for 
our purposes. Median filtering is very effective in remov- 
ing shot noise while preserving the locations of edges. The 
max/min filters are designed to remove morphological fea- 
tures arising from shot noise (see appendix [A")l , 

The filters are defined over neighbourhoods. These 
are often named connectivity or, alternatively, structure 
elements. Image analysis usually deals with regular two- 
dimensional image grids. The most common situation 
for such grids are straightforward 4-connectivities or 8- 
connectivities (see Fig. [3]). When a more arbitrary shape 
is used one usually refers to it as a structure element. 

In the situation of our interest we deal with irregularly 
spaced data, rendering it impossible to use any of the above 
neighbourhoods. It is the Delaunay triangulation which de- 
fines a natural neighbourhood for these situations. For any 
point it consists of its Natural Neighbours, i.e. all points to 
which it is connected via an edge of the Delaunay triangula- 
tion (see Fig. [3j . This may be extended to any higher order 
natural neighbourhood: e.g. a second order neighbourhood 
would include the natural neighbours of the (first order) 
natural neighbours. 

The advantages of following this approach are the same 
as those for the DTFE procedure: the Natural Neighbour fil- 
tering - shortly named NN- median filtering or NN-mm/maa; 



filtering - forms a natural extension to our DTFE based for- 
malism. It shares in the major advantage of being an entirely 
natural and self-adaptive procedure. The smoothing kernel 
is compact in regions of high point concentrations, while it 
is extended in regions of low density. 



3.2.1 Implementation NN Rank- Order filtering 

Implementing the min/max and median Natural Neighbour 
filters within the DTFE method is straightforward. The pro- 
cedure starts with the DTFE density value at each of the 
(original) sample points. These may be the particles in an 
N-body simulation or the galaxies in a redshift survey. For 
each point in the sample the next step consists of the deter- 
mination of the median, maximum or minimum value over 
the set of density values made up by that of the point itself 
and those of its natural neighbours. The new "filtered" den- 
sity values are assigned to the points as the first-order filter 
value. This process is continued for a number of iterative 
steps, each step yielding a higher order filtering step. 

The number of iterative steps of the natural neighbour 
smoothing is dependent on the size of the structure to be 
resolved and the sampling density within its realm. Testing 
has shown that a reasonable order of magnitude estimate is 
the mean number of sample points along the diameter of the 
structure. As an illustration of this criterion one may want 
to consult the low noise and high noise Voronoi models in 
fig. [5] While the void cells of the low noise models contain on 
average 6 points per cell diameter, the void cells of the high 
noise model contain around 16. Fifth-order filtering sufficed 
for the low noise model, 20-th order for the high noise model 
(fig. m and fig.© 

In the final step, following the specified order of the 
filtering process, the filtered density values - determined at 
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Table 1. Parameters of the Voronoi kinematic model realizations: number of cells, cell filling factor, percentages of galaxies within each 
of the morphological components (clusters, filaments, walls, field) and the Gaussian width of clusters, filaments and walls. 



Model 


M 


cell 
filling 
factor 


field 


wall 


(/i- x Mpc) 


filament 


(/i-'Mpc) 


cluster 


Rc 


High noise 


180 


0.500 


50.0 


38.3 


1.0 


10.6 


1.0 


1.1 


0.5 


Low noise 


180 


0.025 


2.5 


16.4 


1.0 


40.6 


1.0 


40.5 


0.5 



the particle positions - are interpolated onto a regular grid 
for practical processing purposes (see sec. 13,1.21 , 

An example of a fifth-order median filtering process is 
shown in the bottom lefthand frame of Fig. [5] The com- 
parison with the original DTFE field (top righthand frame, 
Fig. [2} reveals the adaptive nature of the filtering process, 
suppressing noise in the low-density areas while retaining 
the overall topology of the density field. Figs. [TJa and |8p 
show it in the presence of controlled noise. 



3.3 Markers and False Segment Removal 

Following the NN-median smoothing of the DTFE density 
field, and some minor pixel noise removals, the WVF pro- 
ceeds by identifying the significant minima of the density 
field. These are the Markers for the watershed transform. In 
the case of a cosmological density field the markers are the 
central deep minima in the (smoothed) density field. 

Almost without exception the Markers do not involve 
all minima in a raw unfiltered density field. The minima 
originating from shot noise need to be eliminated. In the 
unfiltered field each regional minimum would correspond to 
a catchment basin, producing over-segmented results: signif- 
icant watershed basins would tend to get subdivided into an 
overabundance of smaller insignificant patches. While most 
of these segments are not relevant a beneficial property of 
the WST is that truely relevant edges constitute a subset of 
the oversegmented segmentation. This notion will be further 
exploited in section [4] 

Once the markers have been selected we compute 
the watershed transform on the basis of an an ordered 
queues algorithm. This process is described in detail in 
|Beucher fc Meverlll993rl . and outlined in appendix [B] The 
process has a few important advantages. It is rather efficient 
because each point is processed only once while it naturally 
involves Watershed by Markers. 



3.4 Hierarchy Merging 

The WVF procedure combines two strategies to remove the 
artefacts generated by Poisson noise resulting from a density 
field discretely sampled by particles or galaxies. 

• the preprocessing of the image such that the insignifi- 
cant minima are removed 

• merging of subdivided cells into larger ones. 



The first strategy involves the previously described recon- 
struction of the density field by DTFE, followed by a com- 
bination of edge preserving median filtering and smoothing 
with the morphological erosion and dilation operators (ap- 
pendix [SJ. In general, as will be argued and demonstrated 
in this study, it provides a good strategy for recovering only 
significant voids. The second strategy involves the merging 
of neighbouring patches via a user-specified scheme. 

Amongst a variety of possibilities we have pursued a 
well known method for merging patches, the watershed hi- 
erarchy. In its original form it assigns to each boundary a 
value dependent on the difference in density values between 
the minima of the neighbouring patches on either side of 
the ridge. We implemented a variant of this scheme where 
the discriminating value is that of the density value inte- 
grated over the boundary. A critical contrast threshold de- 
termines the outcome of the procedure. For an integral den- 
sity value lower than the contrast threshold the two patches 
are merged. If the value is higher the edge is recognized as 
a genuine segment boundary. 

The watershed hierarchy procedure is illustrated in 
Fig. [5ja). An example of its operation is provided by 
Fig. Ob), one of the Voronoi clustering models extensively 
analyzed in the remainder of this study. It depicts the seg- 
mentation resulting from watershed processing of a 5 times 
iteratively NN-median smoothed density field, followed by 
the hierarchical removal of boundaries. The improvement 
compared to the segmentation of a merely 5 times median 
smoothed density field is remarkable (cf. lefthand and right- 
hand panel Fig. [8} . 



3.4.1 Merger Threshold 

In addition to the removal of features on morphological 
grounds, we also have to possibility to remove features on 
the basis of the involved density values. 

In the case of voids we expect that they mature 
as they reach a density defici t of A ~ —0.8 (see e.g 
Shcth fc van de Weygaert |2004| ). Any structures with a 
lower density may be residual features, the diminishing low 
den sity boundaries of th e subvoids which have merged (see 
e.g iDubinski et ~al\ Il993 f ). Various void finding techniques 
do in fact exploit this notion a nd restrict their sea rch to 
regions with A < -0.8 (see e.g. IColberg et al\\20(M ). Note 
that in practice it may also involve noise, of considerable 
significance in these diluted regions. 

A density threshold may indeed be applied within the 
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Figure 6. Frame (a) shows a slice through the original (geometrically defined) Voronoi tessellation. For two different Voronoi clustering 
models defined within this same tessellation, frames (b) and (c) depict the particles within the same slice. Frame (b) shows the low noise 
case with a high density contrast between the voids and walls. Frame (c) is a high noise model with a relatively low contrast between 
voids and walls. 




Figure 7. The density field of the particle distribution in the low noise model (a). Superimposed are the WVF segmentation boundaries. 
The central frame (b) shows the resulting 5-th order median-filtered density field (b). This filtered field is the input for the watershed 
procedure whose segmentation is delineated in frame (c), superimposed on top of the original density field. 




Figure 8. The density field of the particle distribution in the high noise model (a). Superimposed are the WVF segmentation boundaries. 
The central frame (b) shows the resulting 20-th order median-filtered density field. The WVF segmentation of the 5-th order median 
filtered density field, followed by removal of boundaries below a contrast of 0.8, is depicted in frame (c), superimposed on top of the 
original density field. 
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Figure 9. Frame (a): the original (geometric) Voronoi tessellation. Frames (b) and (c): the best recovered WVF segmentation of the 
lownoise (b) and high noise (c) models. 



WVF. This threshold is applied following the watershed 
transform. Any ridges and features with a density contrast 
lower than a specified threshold are removed. The threshold 
A = —0.8 is a natural value of choice. The goal is twofold: 
to suppress noise or spurious features within voids and to 
select out subvoids. 

4 WVF TEST: 

VORONOI CLUSTERING MODEL 

To test and calibrate the Watershed Void Finder 
we hav e applied the WVF t o a K i nematic Voronoi 
Model (Ivan de Weygaert fc Ickd Il989l ; | van de Weygaer tl 
Il99ll . 120021 , 120071 ). In the case of the Voronoi models we 
have exact quantitative information on the location, geom- 
etry and identity of the Voronoi cells, whose interior func- 
tions as the voids in the matter distribution, against which 
we compare the outcome of the WVF analysis. These mod- 
els combine the spatial intricacies of the cosmic web with 
the virtues of a model that has a priori known properties. 
They are particularly suited for studying systematic prop- 
erties of spatial galaxy distributions confined to one or more 
structural elements of nontrivial geometric spatial patterns. 
The Voronoi models offer flexible templates for cellular pat- 
terns, and they are easy to tune towards a particular spatial 
cellular morphology. 

Kinematic Voronoi models belong to the class of 
Voronoi clustering models. These are heuristic models for 
cellular spatial patterns which use the Voronoi tessellation 
as the skeleton of the cosmic matter distribution. The tes- 
sellation defines the structural frame around which matter 
will gradually ass emble during the formation and growth of 
cosmic structure (|Voronoilll90Sl ; IOkabe!l2000h . The interior 
of Voronoi cells correspond to voids and the Voronoi planes 
with sheets of galaxies. The edges delineating the rim of each 
wall are identified with the filaments in the galaxy distribu- 
tion. What is usually denoted as a flattened "supercluster" 
will consist of an assembly of various connecting walls in the 
Voronoi foam, as elongated "superclusters" of "filaments" 
will usually include a few coupled edges. The most outstand- 
ing structural elements are the vertices, corresponding to the 
very dense compact nodes within the cosmic web, rich clus- 
ters of galaxies. We distinguish two different yet complemen- 



tary approaches, Voronoi Element Models and Kinematic 
Voronoi models. The Kinematic Voronoi models are based 
upon the notion that voids play a key organizational role 
in the development of structure and make the Un iverse re- 
semble a soapsud of expanding bubbles llckel (|l984T ). It forms 
an idealized and asymptotic description of the outcome of 
the cosmic structure formation process within gravitational 
instability scenarios with voids forming around a dip in the 
primordial density field. This is translated into a scheme 
for the displacement of initially randomly distributed galax- 
ies within the Voronoi skeleton (see sect IC1I for a detailed 
specification). Within a void, the mean distance between 
galaxies increases uniformly in the course of time. When a 
galaxy tries to enter an adjacent cell, the velocity compo- 
nent perpendicular to the cell wall disappears. Thereafter, 
the galaxy continues to move within the wall, until it tries 
to enter the next cell; it then loses its velocity component 
towards that cell, so that the galaxy continues along a fila- 
ment. Finally, it comes to rest in a node, as soon as it tries 
to enter a fourth neighbouring void. A detailed description 
of the model construction may be found in section [CTI 

To test and calibrate the Watershed Void Finder tech- 
nique we have applied the WVF to a high contrast/low 
noise Voronoi galaxy distribution and a low contrast/high 
noise one. Both concern two stages of the same Kinematic 
Voronoi model, the high noise one to an early timestep with 
a high abundance of field galaxies and the low noise one 
to an advanced stage in which most galaxies have moved 
on towards filament or cluster locations. While the models 
differ substantially in terms of cell filling factor, the under- 
lying geometric pattern remains the same: the position of 
the nodes, edges and walls occupy the same location. Most 
importantly for our purposes: the Voronoi cells, identified 
with the interior of the voids, are the same ones, be it that 
the high noise cells are marked by a substantial population 
of randomly distributed points. 

The model has been set up in a (periodic) box with 141 
h~ Mpc size, and is based on a Voronoi tessellation defined 
by 180 Voronoi cells. In total 128 3 particles were displaced 
following the kinematic Voronoi evolution. Table [1] specifies 
the distinctive parameters defining the model realizations, 
and Fig. [6] shows the particle distribution for the two model 
distributions in a central slice through the model box. 
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4.1 Voronoi Model: Watershed Segmentation 

The density /intensity field is determined by DTFE, yielding 
a 256 3 grid of density values. Fig. [7] contains an example of 
the outcome of the resulting DTFE density interpolation, 
with the contour levels determined according to the descrip- 
tion in section [2] The density map clearly reflects the fila- 
ments and nodes that were seen in the particle distribution. 
The void interiors are dominated by noise, visible as islands 
within a large zero density ocean. 

A direct application of the watershed transform re- 
sults in a starkly oversegmented tessellation (Fig. [7] and 
Fig. |SJ. Amongst the overabundance of mostly artificial, 
noise-related, segments we may also discern real significant 
watersheds. Their boundary ridges (divide lines) are defined 
by filaments, walls and clusters surrounding the voids. Many 
of these genuine voids are divided into small patches. They 
are the result of oversegmentation induced by the noisy Pois- 
son point distribution within the cells. The local minima 
within this background noise will act as individual water- 
shed flood centres marking corresponding, superfluous, wa- 
tershed segments. 

While for a general cosmological distribution it may be 
challenging to separate genuine physical subvoids from ar- 
tificial noise-generated ones, the Voronoi kinematic models 
have the unique advantage of having no intrinsic substruc- 
ture. Any detected substructure has to be artificial, render- 
ing it straightforward to assess the action of the various steps 
intent on removing the noise contributions. 



4-1.1 Smoothing and Segment Merging 

The first step in the removal of insignificant minima consists 
of the application of the iterative natural neighbour median 
filtering process. This procedure, described in sect. 13.21 re- 
moves some of the shot noise in the low density regions. At 
the same time it is edge preserving. The result of five NN- 
median filtering iterations on the high noise version of the 
Voronoi kinematic clustering model is shown in Fig. [7] With 
the exception of a few artificial edges the resulting watershed 
segmentation almost perfectly matches the intrinsic Voronoi 
tessellation. 

Figure [8] shows the result for the high noise version of 
the same Voronoi kinematic clustering model. In this case 
pure NN-median filtering is not sufficient. A much more ac- 
ceptable result is achieved following the application of the 
watershed hierarchy segment merging operation and the re- 
moval of ridges with a density contrast lower than the 0.8 
contrast threshold. 

For both the low-noise and high-noise realizations we 
find that the intrinsic and prominent edges of the Voronoi 
pattern remain in place. Nonetheless, a few shot noise in- 
duced artificial divisions survive the filtering and noise re- 
moval operations. They mark prominent coherent but fully 
artificial features in the noise. Given their rare occurrence we 
accept these oversegmentations as inescapable yet insignifi- 
cant contaminations. 



5 VORONOI CLUSTERING MODEL: 



QUANTITATIVE RESULTS WATERSHED 

The watershed segmentation retrieved by the watershed 
voidfinder is compared with the intrinsic (geometric) 
Voronoi tessellation. The first test assesses the number of 
false and correct WVF detections. A second test concerns 
the volume distribution of the Voronoi cells and the corre- 
sponding Watershed void segments. 

5.1 Datasets 

For our performance study we have three basic models: the 
intrinsic (geometric) Voronoi tessellation, and the low noise 
and high noise Voronoi clustering models (table [TJ. The 
Voronoi clustering models are processed by WVF. In or- 
der to assess the various steps in the WVF procedure the 
models are subjected to different versions of the WVF. 

The second column of Table [2] lists the differently WVF 
processed datasets. These are: 

• Original: the pure DTFE density field, without any 
smoothing or boundary removal, subjected to the watershed 
transform. 

• Mtnmax: only the N~N-min/max filtering is applied to 
the DTFE density field before watershed segmentation. 

• Medn: n iteratations of median natural-neighbour fil- 
tering is applied to the DTFE density field. In all situations 
this includes max/min filtering afterwards. 

• Hierarch: following the watershed transform, on the 
pure non-filtered DTFE density, a density threshold is ap- 
plied. The applied hierarchy threshold level is p/p u = 0.8: 
all segment boundaries with a density lower than S < —0.2 
are removed as physically insignificant. 

• Mednhr. mixed process involving an n times iterated 
median filtered DTFE density field, followed by the water- 
shed transform, after which the segment boundaries below 
the hierarchy threshold 8 < —0.2 are removed. 

Note that the physically natural threshold of A = —0.8 
is not really applicable to the heuristic Voronoi models. On 
the basis of the model specifications the threshold level has 
been set to A = —0.2. 

5.2 Detection Rate 

Each of the resulting segmentations is subjected to a range 
of detection assessments. These are listed in the 3rd to 7th 
column of Table [2] The columns of the table contain respec- 
tively the number of WVF void detections, the amount of 
false splits, the amount of false mergers, the number of cor- 
rectly identified voids, and the correctness measure. While 
the top block contains information on the intrinsic (geomet- 
ric) Voronoi tessellation, the subsequent two blocks contain 
the detection evaluations for the low noise and high noise 
models. 

The false detections are split into two cases. The first 
case we name false splits: a break up of a genuine cell into 
two or more watershed voids. The second class is that of the 
false mergers: the spurious merging of two Voronoi cells into 
one watershed void. The splits, mergers and correct voids 
are computed by comparing the overlap between the vol- 
ume of the Voronoi cell and that of the retrieved watershed 
void. A split is identified if the overlap percentage w.r.t. the 
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Table 2. Quantitative comparison of the original and retrieved voids 



Model 


Parameters 


Voids 


Splits 


Mergers 


Correct 


Correctness 


Intrinsic 




180 












original 


847 












max/min 


259 


82 


3 


118 


66 


Low noise 


med2 


180 


6 


6 


159 


88 




med5 


162 


9 


30 


119 


66 




med20 


136 


20 


80 


33 


18 




original 


4293 












max/min 


3540 













med5 


723 


529 





8 


4 


High noise 


med20 


275 


95 


3 


100 


55 




hierarch 


251 


75 


44 


90 


50 




med5hr 


172 


6 


12 


144 


80 




med20hr 


175 


1 


6 


160 


89 




Figure 10. Lcfthand: the volume distributions for void segments in low-noise models. The histogram shows the intrinsic distribution of 
the Voronoi cell volumes. Superimposed are the inferred volume distribution functions for the WVF segmentations of various Voronoi 
clustering models. The line style of each of the models is indicated in the insert. Righthand: similar plot for a set of noisy Voronoi 
clustering models. 



Voronoi volume is lower than a threshold of 85 percent of the 
overlapping volume. Along the same line, a merger concerns 
an overlap deficiency with respect to the watershed void vol- 
ume. When both measures agree for at least 85 percent a 
void is considered to be correct. The correctness of a certain 
segmentation is the percentage of correctly identified voids 
with respect the 180 intrinsic Voronoi cells. 



5.2.1 Low Noise Model 

Judging by the number of voids in the low noise model, it 
is clear that smoothing or any other selection criterion re- 
main necessary to reduce the number of minima from 850 
to a number close to the intrinsic value 180. The second row 
shows the results for the case when just the maxmin filter 
is applied. This step already reduces the number of insignif- 
icant minima by already 60 percent. It is an indication for 
the local character of the shot noise component. The next 
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Figure 11. Scatter diagram plotting the WVF void segment volumes against the intrinsic geometric Voronoi cell volume. The solid line 
is the linear 1-1 relation. Lefthand: low-noise Voronoi clustering model. Righthand: noisy Voronoi clustering model. 




Figure 12. Scatter diagram plotting the WVF void segment surface area against the intrinsic geometric Voronoi cell volume. The solid 
line is the linear 1-1 relation. Lefthand: low-noise Voronoi clustering model. Righthand: noisy Voronoi clustering model. 



three rows list the results for various iterations of the me- 
dian filtering. With just 2 iterations almost 90 percent of 
the voids are retrieved. Most of the splits are removed at 
2 iterations. This result does not improve with more me- 
dian filtering, even up to 20 iterations this just increases the 
number of mergers as more walls are smoothed away. The 
number of splits also increases as minima begin to merge. 



5.2.2 High noise model 

In general the same conclusion can be drawn for the high 
noise model. Rank-ordered NN-median and NN-min/max 
filters manage to reduce the number of insignificant minima 
by a factor of 80 percent (cf. the number of voids in the 
second and third row). These models attain a correctness of 
approximately fifty percent. Mere rank-ordered filtering is 
evidently insufficient. 
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We also ran a threshold model which did not include 
median filtering. Instead only insignificant boundaries were 
removed. It achieved a recovery of fifty percent. Combining 
both methods (med5hr and med20hr) recovers 80 till 90 per- 
cent of the voids. The succes rate may be understood by the 
complementarity of both methods: while the median filter- 
ing recovers the coherent structures, the thresholding will 
remove those coherent walls that are far underdense. 

The translation to a cosmological density field is 
straightforward. The rank-ordered filtering ensures that in- 
significant minima are removed and that the watershed will 
pick up only coherent boundaries. Thresholding is able to 
order these walls by significance and to remove the very un- 
derdense and insignificant walls. 

5.3 Volume Comparison 

In Fig. [10] we compare the distribution of the void volumes. 
The histogram shows the distribution of equivalent radii for 
the segment cells, 

r - {¥■ w 

The solid line histogram shows the (geometric) volume dis- 
tribution for the intrinsic Voronoi tessellations. On top of 
this we show the distributions for the various (parameter- 
ized) watershed segmentation models listed in Table [2] Not 
surprisingly the best segmentations have nearly equivalent 
volume-distributions. For the lownoise models this is med2 
(lefthand), for the highnoise models med20hr (righthand). 
This conclusion is in line with the detection rates listed in 
Tabled 

The visual comparison of the intrinsic geometric 
Voronoi tessellations and the two best segmentations - med2 
for the lownoise model and med20hr for the highnoise ver- 
sion - confirms that also the visual impression between these 
watershed renderings and the original Voronoi model is very 
much alike. 

We have also assessed the cell-by-cell correspondence 
between the watershed segmentations and the Voronoi 
model. Identifying each watershed segment with its origi- 
nal Voronoi cell we have plotted the volume of all watershed 
cells against the corresponding Voronoi cell volumes. The 
scatter plots in Fig. [TTJ form a convincing confirmation of 
the almost perfect one-to-one relation between the volumes 
derived by the WVF procedure and the original volumes. 
The only deviations concern a few outliers. These are the 
hierarchy merger segments for which the watershed volumes 
are too large, resulting in a displacement to the right. 

5.4 Surface Comparison 

While the volumes occupied by the watershed segments in 
Fig. [9] do overlap almost perfectly with that of the original 
Voronoi cells, their surfaces have a more noisy and erratic 
appearance. This is mostly a consequence of the shot noise 
in the (DTFE) density field, induced by the noise in the 
underlying point process. The crests in the density field are 
highly sensitive to any noise, 

In addition to assess the impact of the noise on the sur- 
faces of the watershed segments we compared the watershed 



segement surface areas with the Voronoi cell surface areas. 
The results are shown in Fig. 1121 We tested the lownoise 
med2 and the highnoise med20hr. In both cases we find a 
linear relationship between the watershed surface and the 
genuine Voronoi surface area. Both cases involve consider- 
ably more scatter than that for the volumes of the cells. In 
addition to an increased level of scatter we also find a small 
be it significant offset from the 1-1 relation. The slope of 
the lownoise model is only slightly less than unity, the high- 
noise model slope deviates considerably more. These offsets 
do reflect the systematically larger surface areas of the wa- 
tershed segments, a manifestation of their irregular surfaces. 
It is evident that the level of irregularity is more substantial 
for the highnoise than for the lownoise reconstructions (cf. 
Fig.©. 

The scatter plots do also reveal several cells with huge 
deviations in surface area. Unlike expected there is no sys- 
tematic trend for smaller cells to show larger deviations. 
Some of the small deviating cells can be recognized in Fig. [9] 
as highly irregular patches. The large deviant cells corre- 
spond to watershed segments which as a result of noisy 
boundaries got wrongly merged. 

While the irregularity of the surface areas forms a good 
illustration of the noise characteristics of the watershed 
patches, for the purpose of void identification it does not 
pose a serious problem. Smoother contours may always be 
obtained by applying the flooding process on a properly 
smoothed field. Some suggestions for this may be achieved 
follows in the discussion. 



6 DISCUSSION AND PROSPECTS 

The WVF void finder technique is based on the watershed 
transform known from the field of image processing. Voids 
are identified with the basins of the cosmological mass dis- 
tribution, the filaments and walls of the cosmic web with the 
ridges separating the voids from each other. Stemming from 
the field of mathematical morphology, watershed cells are 
considered as patches locally minimizing the "topographic 
distance" . 

The WVF operates on a continuous density field. Be- 
cause the cosmological matter distribution is nearly always 
sampled by a discrete distribution of galaxies or N-body par- 
ticles, the first step of the WVF is to translate this into a 
density field by means of the Delaunay Tessellation Field 
Estimator (DTFE). Because the WVF involves an intrinsi- 
cally morphological and topological operation, the capability 
of DTFE to retain the shape and morphology of the cosmic 
web is essential. It guarantees that within this cosmological 
application the intrinsic property of the watershed trans- 
form to act independent of scale, shape, and structure of a 
segment is retained. As a result, voids of any scale, shape 
and structure may be detected by WVF. 

In addition to the regular watershed transform the 
WVF needs to invoke various operations to suppress (dis- 
creteness) sampling noise. In addition, we extend the water- 
shed formalism such that the WVF will be capable of analyz- 
ing the hierarchy of voids in a matter distribution, i.e. iden- 
tify how and which smal l scale voids are embedded within 
a voi d on larger scales (|Platen. van de Weygaert fc Jones! 
2007). Markers indicating significant void centers, false seg- 
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ment removal by Hierarchy Merging and Natural Neighbour 
filtering all affect an efficient noise removal. Hierarchy Merg- 
ing manages to eliminate boundaries between subvoids. Nat- 
ural Neighbour median filtering, for various orders, is an es- 
sential new ingredient for highlighting the hierarchical em- 
bedding of the void population. It allows a natural selection 
of voids, unbiased with respect to the scale and shape of 
these structures. The voids that persist over a range of scales 
are expected to relate to the voids that presently dominate 
the cosmic matter distri bution. In other words, WVF pr e- 
serves the void hierarchy ISheth fc van de Weygaertl {2004). 

The present work includes a meticulous qualitative and 
quantitative assessment of the watershed transform on the 
basis of a set of Voronoi kinematic models. These heuristic 
models of spatial weblike or cellular galaxy or particle dis- 
tributions facilitate the comparison between the void detec- 
tion of the WVF and that of the characteristics of the cells 
in the original and intrinsically known Voronoi tessellation. 
It is found that WVF is not only succesfull in reproducing 
the qualititative cellular appearance of the Voronoi models 
but also in reproducing quantitative aspects like an almost 
perfect 1-1 match of cell size with WVF segment volume 
and the corresponding void size distribution. 

We foresee various possible improvements of the WVF. 
These concern in particular the identification of signifi- 
cant edges. One po ssibility is that extension proposed by 
iNguven et all ([2003:) , in which not only the "topographic 
costs" but also the lengths of the contours should be min- 
imized. The length minimization will result in smoother 
boundaries. Additional improvements may be found in 
better filtering procedures in order to facilitate studies 
of hierarchically structured patterns. We expect consid- 
erable improvements by anisotropic diffusion techniques 
(|Black fc Marimontl I1998T ) and are currently implementing 
these in the WVF code. 

Given the results of our study, we are confident for ap- 
plying WVF to more elaborate and realistic simulations of 
cosmic structure formation and on large galaxy redshift sur- 
veys. The analysis of a set of GIF cosmological simulations 
will be presented in an upcoming paper. 
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APPENDIX A: 

MATHEMATICAL MORPHOLOGY 

Appendix [X] and [B] provide some formal concepts and no- 
tations necessary for appreciating the watershed transform. 

The WVF formalism introduced in this study is largely 
based upon concepts and techniques stemming from the field 
of image analysis. Although they are used within the con- 
text of the morphological analysis of spatial patterns in cos- 
mological density fields the presentation is in terms of the 
original image analysis nomenclature. In this we remind the 
cosmology reader to translate "image" into "density field", 
"basin" into "void interior", etc. 

Mathematical Morphology (MM) is the field of image 
analysis which aims at characterizing an image on the ba- 
sis of its geometrical structure. It provides techniques for 
extracting image components which are useful for repre- 
sentation and description and was originally developed by 
G. Mathero n and J. Serra. For more details we r efer t o 
Serral dl983l), also see iHeiimansI l| 19941 ): iMatheronl l|l975h : 
Meyer fe Beucherl (|l990r7 7 It involves a set-theoretic method 
of image analysis and incorporates concepts from algebra 
(set theory, lattice theory) as well as from geometry (trans- 
lation, distance, convexity). Applications of mathematical 
morphology may be found in a large variety of scientific 
disciplines, including material science, medical imaging and 
pattern recognition. 



Al Images 

A cosmological density field /(x) may be mapped onto an 
image T, 

/(x) -» ^(x). (Al) 

The image T is a function in n-dimensional lattice space 
(usually n — 3 or n = 2) Z™, 

T : Z™ -f Z . (A2) 

Although in principle images may be continuous, in practice 
they usually attain a finite number of discrete values. Two 
important classes are: 
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Figure Al. Illustration of the effects of a few essential operators on a binary image (top column) and a grayscale image (bottom 
row). The original image is the one in the lefthand frame. The central frame contains the image following a opening operation, while 
the righthand frame shows the effect of closing. The circle at the lower lefthand corner of the originals represent the circular structure 
element . 



• Binary image: 
image with only 2 intensity values (Fig. lAll top lefthand), 



(A3) 



We follow the convention to identify the binary image by 
the set X C Z n for which T = 1. 

• Grayscale image: 
image with a discrete number N of values (Fig. IA1I bottom 
lefthand), 



T = i = l,...,N} 



(A4) 



Mathematical morphology was originally developed for bi- 
nary images, and was later extended to grayscale images. 



A2 Erosion and Dilation 

The two basic operators of Mathematical morphology are 
Erosion and Dilation of a binary image X. In order to de- 
fine these operators we need to invoke the translation and 



reflection of a set B, 

translation : B z 
reflection : B 



{y I V = b + z 
{y\y = -b 



V b € B} 

V b 6 B} 



(A5) 

The dilation or erosion of a binary image X by a structuring 
element B identifies whether the translated set B z has an 
overlap with or is contained in a certain part of X, 



dilation 



X®B 
X Q B 



B z n x / 1} 

B z C X} 



(A6) 



In other words, dilation consists of the Minkowski addition 
(©) of the binary image X with a structuring element B 
while erosion is the Minkowski substraction (©) with B. A 
structuring element may be any object in Z n . An example 
is the circle which functioned as a structuring element in 
Fig. IA1I Erosion and dilation have a number of properties: 



Translation invariance 
Global scaling invariance 
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• Addition: 

if X CY -> X ® ACY ® A (A7) 

• Complementarity: 

erosion of a set is dilation of complement 
(and vice versa) 

• Adjunction relationship: 

Y © A C X <^ Y CXQA. (A8) 

The complementarity and adjunction relationship are two 
aspects of the existing duality between erosion and dilation. 
In general erosion and dilation induce a loss of information. 
Erosion followed by dilation, or vice versa, will only result 
in a restoration of the original image if the sets X, Y and A 
are convex. In fact, various combinations of the erosion and 
dilation operators do result in new operators. 



of the binary image definitions (eqn. IA5)) implies the follow- 
ing definition wrt. a grayscale image, 

T®B = s\xp{F(x + b), x£X,beB} s 
TQB = m{{T(x + b), x€X,beB} 1 ' 

The effect of erosion on a grayscale image is the shrinking 
of the bright regions. Bright spots smaller than the struc- 
turing element B disappear completely while valleys (dark) 
expand. Dilation has the opposite effect: dark regions shrink 
while bright regions expand. It illustrates the duality be- 
tween erosion and dilation. 

For our purposes, this formal definition translates into 
the following practical implementation for 2-D grayscale 
images. Given a grayscale image T(A) with grid elements 
a(i,j), 



A3 Opening and Closing 

The most straightforward combination of dilation and ero- 
sion is that of the consecutive application of an erosion and 
a dilation. This introduces two new operators, the Opening 
and Closing operators. On a binary image they have the ef- 
fects shown in Fig. IA1I (top centre, top righthand). Opening 
amends caps, removes small islands and opens isthmuses. 
Closing, on the other hand, closes channels, fills small lakes 
and (partly) the gulfs. An additional combination of erosion 
and dilation is subtraction of the first from the latter, called 
a morphological gradient. 

Formally, an opening is an erosion followed by a dila- 
tion, a closing a dilation followed by erosion. 

Opening: A B = [(X 6 B) © B] 

Closing: $ s = [{X (B B) Q B] ( ' 

Characteristics of the opening and closing operators are: 

• Increasing 

• Idempotent: 

applying the operator twice yields the same output 

• Opening is anti-extensive: 

Apr) C X (A10) 

• Closing is extensive: 

X C $(A) (All) 

Note that the extensivity and/or anti-extensivity of opera- 
tors define the prime conditions for a morphological operator. 



A4 Grayscale Images 

The morphological operators which we discussed above can 
be generalized to grayscale images. 

A grayscale image is composed of subsets Si(J-), 

Si(F) = {x\x G Z" : T{x) i}. (A12) 

with Si+^T) C Si(T) and S\ (T) the support of the full im- 
age. The erosion and dilation of a grayscale image involves 
their application to each individual subset Si(J-). Extension 



T®B = maX(ij) eB {a[r + i, s + j] + b\i,j]} 

V[r, s] e A 

TQB = min (M)es {a[r + i,s + j] + b[i,j]} 

V[r, s] e A 
(A14) 

As in the case of binary images new operators may be 
defined through combinations of erosions and dilations. The 
closing and opening operators are defined in exactly the 
same way as that for the binary images. Their effect is shown 
in the lower column of Fig. IA1I The morphological gradient 

g = (T © B) © (FeB) (A15) 

is a dilation minus erosion operation. The gradient operator 
is often used in object detection because an object is usu- 
ally associated with a change in grayscale with respect to 
the background. A variety of additional operators involving 
openings and closings may be defined. Interesting ones are 
the granulometries, a sequence of erosions with increasing 
scale, and distance transforms. 



APPENDIX B: 
WATERSHED TRANSFORM 

The segmentation of images is defined on the basis of a dis- 
tance criterion, referring to the concept of distance between 
subsections of an image. 

Bl Distance 

For appreciating the concept of watershed segmentation we 
consider two distance concepts, the geodesic and topographic 
distance. Geodesic distances are used in the case of binary 
images while topographic distances form the basis for the 
segmentation of grayscale images. 

Let X C Z™ be a set and x and y two points in an 
n-dimensional lattice space Z™. We may define 

Bl.l Geodesic Distance 

The geodesic distance dx(x,y) is the length of the short- 
est (geometric) path in X connecting x and y (see lefthand 
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frame, Fig. IB1|I . Accordingly, the distance between two sub- 
sets A and B in X may be defined as follows. Considering 
the set of all paths between any of the elements of A and 
those of B, the distance dx(A,B) between A and B is de- 
fined to be the minimum length of any of these paths. 

Based upon the concept of dx(A, B) one may formulate 
a distance function of a set Y G Z™. For each point y £ 
Y, the distance d(y,Y) to its complement Y is computed. 
The distance function T>(y) is the resulting map of distance 
values dx(y, Y). Regions whose distance dx(y, Y) is at least 
Ti can be identified and equated by erosion of Y with a disk 
of radius Ti (defining a set IZi). Each of these regions forms 
a section St>4 (see sec. IA4|) . 

&D,t = V{y)®Ki, (Bl) 

in which IZi is the disk of radius rt. The map T>(y) may be 
regarded as a stack of these sections. For illustration the dis- 
tance transform of the the binary image Fig. lAll is depicted 
in the central frame of Fig IB 1 1 

B1.2 Topographic Distance 

The topographic distance T(x, y) between two points in Z™ 
is defined with respect to the image map T. Taking the limit 
of a continuous map T , the topographic distance from x to 
y is the path which attains the minimum length through the 
"image landscape", 

T{x,y) = inf y \VF(j(s))\ds. (B2) 

In this definition, the integral denotes the image pathlength 
^-"(7) along all paths y(s) in the set of all possible paths, V. 
This concept of distance is related to the geodesies of the 
surface T: the path of steepest descent, specifing the track 
a droplet of water would follow as it flows down a mountain 
surface. 



B2.2 Skeleton 

The boundary set in X consists of those points which do 
belong to X yet are not contained in any of the zones of 
influence Zp(Yi). These boundary points define the geodesic 
Skeleton K. of T, 

Kr = X\Z T . (B6) 

In Fig. IB 1 1 frighthand frame) the skeleton in set X is outlined 
by white lines. The skeleton is superimposed on its defining 
distance function landscape, its values indicated by a red 
colour gradient scheme (the corresponding landscape profile 
is depicted in the central frame). We should point out that 
here we follow the definition of mathematical morphology al- 
though the name of skeleton of the large scale cosmic matter 
distribu tion has been used for different be it related c oncept s 
(see e.g. Ivan de Weygaertlll99ll ; iNovikov ei"ai1l2006l '). 

It is interesting to note that if we restrict the subsets Yi 
to single points the skeleton naturally evolves into a (first- 
order) Voronoi Tessellation. It is the definition of the con- 
cept of skeleton K, within the specific context of grayscale 
images which brings us to the definition of the Watershed 
Transform (see next section). 

B3 The watershed transform: Algorithms 

Grayscale images consist of a finite number of discrete lev- 
els. This results in a slightly more complicated situation for 
its segmentation. In the case of a binary image an image is 
segmented on the basis of a geodesic distance. For the seg- 
mentation of grayscale images the distance concept needs to 
be generalized to that of topographic distance. 

Definition: Each watershed basin is the collection of points 
which are closer in topographic distance to the defining min- 
imum then to any other minimum. 



B2 Segmentation 

Based on the specific definition of distance, we may segment 
a binary image T through the identification of the zones of 
influence of well-defined subsets of T. In general the binary 
image T contains n connected subsets Yi [i = 1, . . . , n) (the 
black regions of Fig. IA1|) and a set X (white region) which 
contains all points x which do not belong to any of the sub- 
sets Yi, i.e. 

X = \x&F\x i \JyA . (B3) 



B2. 1 Zone of Influence 

The geodesic zone of influence Zp(Yi) of a subset Yi G T is 
the set of all the points x £ X that are stricktly closer to Yi 
than to any other subset Yj, j 7^ i, 

ZrQTi) = {x€X\d x {x,Y) <dx{x,Y 3 ), V j + i } (B4) 

The Zone of Influence Z of T is the union of all influence 
zones of Z^{Yi), 

Z T = |J ZrQQ (B5) 



The literature is replete with algorithms for the con- 
struction of the watershed transform of an image into its 
constituting watershed segments. They may be divided into 
two classes. One class simulates the watershed basin immer- 
sion process. The second aims at detecting the watershed 
skeleton on the basis of the distance. 

B3.1 Watershed by Immersion 

Watershed by immers ion was introduced and defined by 
IVincent fc Soillel (|l99ll ). The first step of the procedure con- 
cerns the identification of the minima. Formally, a minimum 
is a plateau at altitude h from which it is impossible to reach 
a point of lower height. Starting with the lowest grayscale 
level hmin and recursively proceeding to the highest level 
h m ax the algorithm allocates the zone of influence of each 
minimum by gradually filling up the surrounding catchment 
basin. 

At a particular grayscale level i, with altitude hi, the 
algorithm has hypothetically inundated a landscape region 
with an altitude .F(x) ^ hi. The total of inundated area, 

Sr{T) = {x e TT\T(x) < i}, (B7) 

is the complement of the section Si. Having arrived at level 



© 2007 RAS, MNRAS 000. 011241 



20 Platen, van de Weygaert & Jones 




Figure Bl. These three images illustrate the concept of distance and segmentation. Left: the geodesic distance between points x and 
y is depicted within a binary image. Centre/Right: a landscape defined by a distance function (centre) and the resulting segmentation 
and skeleton (white line, righthand frame). 



i and proceeding to level i + 1 the algorithm has three pos- 
sibilities, 

1 encounter a new minimum (at level i + 1) 

2 add new points to existing catchment basins (condition: 
points connected to only one existing basin). 

3 encounter new points that belong to several basins 

Situation (1) signals the event that a new minimum becomes 
active in the image. The second option concerns the exten- 
sion of an existing basin by an additional collection of points. 
These are points identified at level (i + 1) which find them- 
selves within the realm of a single basin existent at level i. 
They belong to the zone of influence of Si and find them- 
selves embedded in its extended counterpart Si+i at level 
(i + 1). In situation (3) more than one basin may be con- 
nected to the stack Si+i. Its correct subdivision is deter- 
mined by computing the influence zones of all connected 
basins. Defining Bi(T) to be the union of all catchment 
basins at level i the union of catchment basins at level (i+ 1) 
becomes 



B i+1 {H) 



u Bi 



(B8) 



The watershed procedure may be viewed as iteratively com- 
puting the zone of influence at each new grey scale level. 

Following the rationale above the final "immersion" def- 
inition for the watershed W of the image T within a domain 
X is 



W(H = X \ B hn 



(B9) 



On completion of the procedure the union of points attached 
to every minimum m in X is equal to the union of catch- 
ment basins, Bh ma:c ■ The skeleton remains as the watershed 
segmentation. 



B3.2 Watershed by Topographic Distance 

The alternative strategy for determining the water- 
shed transform is that of following the strict defini- 
tion of segmentation by m inimum topographic distance. 
iRoerdink fc Meiisterl ((2000) give a summary of the most 
notable schemes. These algorithms seek to find all points 



(pixels) whose topographic distance to a particular marker - 
ie. a significant minimum in the density field - is the shortest 
amongst that to all other markers in the image. 

The formalism bears some resemblance to Dijkstra's 
graph theoretical problem of tracing the shortest path forest 
in a point distribution. Based on this similarity an image is 
seen as a connected (di)graph in which the pixels of the im- 
age function as the nodes of the graph. Each point p is reach- 
able from each other point p' via the graph's edges. The lat- 
ter usually define a network on the basis of 4— connectivities 
or 8-connectivities. 

The shortest path between two points (nodes) p and p' 
is found by traversing the graph and keeping track of the 
walking cost. Critical for the procedure is the assignment of 
a proper measure of cost to each path. By definition it should 
be a non-negative increasing function and be related to the 
definition of topographic distance (eq. IB2|I . This suggests 
the use of the quantity 



H(p,p) — max 



[ Hp) -Hp') 

I d(p,p') 



(BIO) 



the maximum slope linking the two pixels p and p . This 
leads to the following cost function for the link between two 
neighbouring pixel p and pixel p , 



H(p,p')d(p, P ') 



Hp) > Hp') 



c(p,p') = I T'{p',p)d{p, P ') Hp) < Hp') 

H(p,p')+H(p', P ) d{ppl) m=T{p>) 

(Bll) 

The total cost C for a path "f(pi,P2, ■■■Pn) connecting any 
two points pi and p„ via the points {p2, . . . ,p n -i} will then 
simply be the sum 



C(pi,p n ) = ^C(p,,p l _i). 



(B12) 



The topographic distance T(pi,p n ) is the infimum of 
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C J (pi,p n ) over all paths connecting pi and p n . 

T(p uPn ) = inf C 7 0i,p„) (B13) 

Given this definition for the topographic distance within a 
grayscale image we can pursue the segmentation process de- 
scribed in section [B2l ultimately yielding the watershed seg- 
mentation. 



B4 Ordered Queues Algorithm 

We follow the water shed transform algorithm by 
iBeucher fc Meverl (|l993l ). Their method implicitly in- 
corporates the concept of Markers. These markers are 
the minima used as sources of the watershed flooding 
procedure. As such they form a select subgroup amongst all 
minima of an image T . 

The code for the watershed proceudre involves the fol- 
lowing steps: 

• Initialization All pixels of the cube are initialized 
and tagged to indicate they have not yet been processed. 
Each grayscale level is allocated a queue and all pixels are 
attached to the queue corresponding to their level. 

• Minima Each minimum plateau is tagged by a unique 
"minimum tag" . The pixels corresponding to a minimum are 
inserted into the corresponding queue. 

• Flooding All pixels in the grayscale level queues are 
processed, starting at the lowest grayscale level. Unless a 
pixel is surrounded by a complex of unprocessed neighbours 
it will be assigned to the queue of the corresponding mini- 
mum. Pixels which also border another minimum obtain a 
boundary tag. 

• Final Stage For any grayscale level the flooding stops 
when the queue has emptied. The procedure continues with 
processing the pixels in the queue for the next grayscale 
level. The process is finished once all level queues have been 
emptied. 



APPENDIX C: 

KINEMATIC VORONOI MODELS 

Voronoi Clustering Models are a class of heuristic models 
for ce llular dist ributions o f matter Ivan dc Wcyga ert fc 
(1989); van do WeveaertJ <|l99ll 120021 . l2007h . They use the 
Voronoi tessellation as the skeleton of the cosmic matter 
distribution, identifying the structural frame around which 
matter will gradually assemble during the emergence of cos- 
mic structure. The interior of Voronoi cells correspond to 
voids and the Voronoi planes with sheets of galaxies. The 
edges delineating the rim of each wall are identified with the 
filaments in the galaxy distribution. The most outstanding 
structural elements are the vertices, corresponding to the 
very dense compact nodes within the cosmic web, the rich 
clusters of galaxies. 

We distinguish two different yet complementary ap- 
proaches. One is the fully heuristic approach of Voronoi El- 
ement models. They are particularly apt for studying sys- 
tematic properties of spatial galaxy distributions confined to 
one or more structural elements of nontrivial geometric spa- 
tial patterns. The second, supplementary, approach is that 
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Figure CI. Schematic illustration of the Voronoi kinematic 
model. Courtesy: Miguel Aragon. 



of the Voronoi Kinematic models, which attempts to "simu- 
late" foamlike galaxy distributions on the basis of simplified 
models of the evolution of the megaparsec scale distribution. 

The Voronoi Kinematic Model is based upon the notion 
that voids play a key organizational role in the development 
of structure and m ake the Un iverse resemble a soapsud of 
expanding bubbles llckei |l984) . It forms an idealized and 
asymptotic description of the outcome of the cosmic struc- 
ture formation process within gravitational instability sce- 
narios with voids forming around a dip in the primordial 
density field. For plausible structure formation scenarios, 
most notably the concordance ACDM cosmology, this evo- 
lution will proceed hierarchically . A detailed assessment of 
the r esulting void hierarchy by l|Sheth fc van de Weygaert] 
|2004) demonstrated that this leads to a selfsimilarly evolv- 
ing peaked void size distribution. By implication, most voids 
have comparable sizes and excess expansion rates. The ge- 
ometrically interesting implication is that the asymptotic 
limit of the "peaked" void distribution degenerating into one 
of only one characteristic void size. It yields a cosmic matter 
distribution consisting of equally sized and expanding spher- 
ical voids, a geometrical configuration which is precisely that 
of a Voronoi Tessellation. This is translated into a scheme 
for the displacement of initially randomly distributed galax- 
ies within the Voronoi skeleton (see sect. ICf 1 for a detailed 
specification). Within a void, the mean distance between 
galaxies increases uniformly in the course of time. When a 
galaxy tries to enter an adjacent cell, the velocity compo- 
nent perpendicular to the cell wall disappears. Thereafter, 
the galaxy continues to move within the wall, until it tries 
to enter the next cell; it then loses its velocity component 
towards that cell, so that the galaxy continues along a fila- 
ment. Finally, it comes to rest in a node, as soon as it tries 
to enter a fourth neighbouring void. 
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Cl Initial Conditions 

The initial conditions for the Voronoi galaxy distribution 
are: 

• Distribution of M nuclei, expansion centres, within the 
simulation volume V. The location of nucleus Tf I IS ym - 

• Generate N model galaxies whose initial locations, x„o 
(n = 1,...,N), are randomly distributed throughout the 
sample volume V. 

• Of each model galaxy n determine the Voronoi cell V a 
in which it is located, ie. determine the closest nucleus j a . 

All different Voronoi models are based upon the displace- 
ment of a sample of N "model galaxies" . The initial spatial 
distribution of these N galaxies within the sample volume V 
is purely random, their initial locations x„o (n — 1, . . . , N) 
defined by a homogeneous Poisson process. A set of M nu- 
clei within the volume V corresponds to the cell centres, or 
expansion centres driving the evolving matter distribution. 
The nuclei have locations y m (m = 1, . . . , M). 

Following the specification of the initial positions of all 
galaxies, the second stage of the procedure consists of the 
calculation of the complete Voronoi track for each galaxy 
n — 1, . . . ,N (sec. IC2[l . Once the Voronoi track has been 
determined, for any cosmic epoch t one may determine the 
displacement x n that each galaxy has traversed along its 
path in the Voronoi tessellation (sec. IC2|l . 

C2 Voronoi Tracks 

The first step of the formalism is the determination for each 
galaxy n the Voronoi cell V a in which it is initially located, 
ie. finding the nucleus j a which is closest to the galaxies' 
initial position x n o- 

In the second step the galaxy n is moved from its ini- 
tial position x„o along the radial path emanating form its 
expansion centre j a , ie. along the direction defined by the 
unity vector e na . Dependent on how far the galaxy is moved 
away from its initial location x„o - set by the radius of ex- 
pansion R n to be specified later - the galaxies' path x n (see 
Fig. ICip may be codified as 

Xn — yet S ncK -(- S nQ ^ -)- Snafiy 

(Cl) 

— ya "T S na G na -\- S na fi(i na {3 "T S na p^G na f3y 

in which the four different components are: 

• e na : unity vector path within Voronoi cell V Q 

• &naf)- unity vector path within Voronoi wall E Q/ g 

• e na/ 3 7 : unity vector path along Voronoi edge A aj 9 7 

• Vertex H a/ 3 7( s 

The identity of the neighbouring nuclei j a , jp, j 7 and j$, 
and therefore the identity of the cell V a , the wall £ a /3, the 
edge A a/ 3 7 and the vertex H a/ 3 7 ,5, depends on the initial loca- 
tion x„o of the galaxy, the position y a of its closest nucleus 
and the definition of the galaxies' path within the Voronoi 
skeleton. 

The cosmic matter distribution at a particular cos- 
mic epoch is obtained by calculating the individual dis- 
placement factors {snait), Snaftit), s na fs-y{i)) for each model 



galaxy. These are to be derived from the global "void" ex- 
pansion factor R(t). This factor parameterizes the cosmic 
epoch and specifies the (virtual) radial path of the galaxy 
from its expansion centre j a . 

At first, while still within the cell's interior, the galaxy 
proceeds according to 

s na {t) = x n (t) - y a = R(t) |x n0 -y a \e na . (C2) 

As a result within a void the mean distance between galaxies 
increases uniformly in the course of time. Once the galaxy 
tries to enter an adjacent cell jp and reaches a Voronoi wall, 
i.e. when R(t) |x„o — y a > v n , the galaxy's motion will be 
constrained to the radial path's component within the wall 
E Q( 3. The galaxy moves along the wall until the displace- 
ment supersedes the extent of the path within the wall and 
it tries to enter a third cell j 7 , i.e. when s na/ s(t) > a n . Sub- 
sequently, it moves along A a /3 7 until it comes to rest at the 
node H a/ 3 7 a as soon as it tries to enter a fourth neighbouring 
void js when s nQ( 3 7 > A„. 

A finite thickness is assigned to all Voronoi structural 
elements. The walls, filaments and vertices are assumed to 
have a Gaussian radial density distribution specified by the 
widths i?w of the walls, Rf of the filaments and Rv of the 
vertices. Voronoi wall galaxies are displaced according to 
the specified Gaussian density profile in the direction per- 
pendicular to their wall. A similar procedure is followed for 
the Voronoi filament galaxies and the Voronoi vertex galax- 
ies. As a result the vertices stand out as three-dimensional 
Gaussian peaks. 



APPENDIX D: 

DTFE RECONSTRUCTION PROCEDURE 

For a detailed spe cification of th e DTFE density field pro- 
cedure we refer to ISchaad (|2007h . In summary, the DTFE 
procedure for density field reconstruction from a discrete set 
of points consists of the following steps: 

• Point sample 

Given that the point sample is supposed to represent an 
unbiased reflection of the underlying density field, it needs 
to be a general Poisson process of the (supposed) underlying 
density field. 

• Boundary Conditions 

The boundary conditions will determine the Delaunay and 
Voronoi cells that overlap the boundary of the sample vol- 
ume. Dependent on the sample at hand, a variety of options 
exists: 

+ Empty boundary conditions: 
outside the sample volume there are no points. 

+ Periodic boundary conditions: 
the point sample is supposed to be repeated periodically 
in boundary boxes, defining a toroidal topology for the 
sample volume. 

+ Buffered boundary conditions: 
the sample volume box is surrounded by a bufferzone filled 
with a synthetic point sample. 

• Delaunay Tessellation 

Construction of the Delaunay tessellation from the point 
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Figure C2. Evolution of galaxy distribution in the Voronoi kinematic model. A sequel of 3 consecutive timesteps within the kinematic 
Voronoi cell formation process, proceeding from left to right, and from top to bottom. The depicted boxes have a size of 100/i _1 Mpc. 
Within these cubic volumes some 64 Voronoi cells with a typical size of 25fo _1 Mpc delineate the cosmic framework around which some 
32000 galaxies have aggregated. Taken from a total (periodic) cubic "simulation" volume of 200/i — ^^Mpc containing 268,235 "galaxies". 



sample. While we also s t ill use the Voronoi-Delaunay code of 
( van de Wevgaert]|l99ll.ll994r) , at present there is a number 
of efficient library routines available. Particularly notewor- 
thy is the CGAL initiative, a large library of computational 
geometry routines 1 

• Field values point sample 

The estimate of the density at each sample point is the 
normalized inverse of the volume of its contiguous Voronoi 
cell Wj of each point i. The contiguous Voronoi cell of a 
point i is the union of all Delaunay tetrahedra of which 
point i forms one of the four vertices. We recognize two 
applicable situations: 

- uniform sampling process: 
the point sample is an unbiased sample of the underlying 
density field. Typical example is that of TV-body simulation 
particles. For D-dimensional space the density estimate is, 



j3(x0 = (! + £>) 



v{m) 



(Dl) 



with Wi the weight of sample point i, usually we assume 
the same "mass" for each point. 

- systematic non-uniform sampling process: 
sampling density according to specified selection process. 
The non-uniform sampling process is quantified by an a pri- 
ori known selection function V( x )- This situation is typi- 
cal for galaxy surveys, ^>(x) may encapsulate differences in 
sampling density ip(a,5) as function of sky position (a, S), 
as well as the radial redshift selection function ip(r) for 
magnitude- or flux-limited surveys. For D-dimensional space 
the density estimate is , 



p( Xi ) = (1 + D) 



(D2) 



1 CGAL is a C++ library of algorithms and data structures for Com- 
putational Geometry, see www . cgal . org 



• Field Gradient 

Calculation of the field gradient estimate V/| m in each 
D-dimensional Delaunay simplex m (D = 3: tetrahedron; 
D = 2: triangle) by solving the set of linear equations 
for the field values fi at the positions of the (D + 1) 
tetrahedron vertices, 



v/|. 



fo fi h h 



. r ri r 2 r 3 



(D3) 



Evidently, linear interpolation for a field / is only meaningful 
when the field does not fluctuate strongly. 

• Interpolation. 

The final basic step of the DTFE procedure is the field in- 
terpolation. The processing and postprocessing steps involve 
numerous interpolation calculations, for each of the involved 
locations x. Given a location x, the Delaunay tetrahedron 
m in which it is embedded is determined. On the basis of the 
field gradient V/| m the field value is computed by (linear) 
interpolation, 



/(x) = /(xO + V/L -(x-Xj 



(D4) 



In principle, higher-order interpolation procedures are also 
possible. Two relevant procedures are: 

- Spline Interpolation 

- Natural Neighbour Interpolation 

For NN-interpo la tion see IWatsonl (1 19921): 
iBraun fc Sambridgel l|l995l ): ISukumarl (199^ and lOkabel 
1 20001 ). Implementation of Natural neighbour interpolations 
is presently in progress. 

• Processing. 

Though basically of the same character, for practical pur- 
poses we make a distinction between straightforward pro- 
cessing steps concerning the production of images and simple 
smoothing filtering operations and more complex postpro- 
cessing. The latter are treated in the next item. Basic to the 
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processing steps is the determination of field values following 
the interpolation procedure(s) outlined above. Straightfor- 
ward "first line" field operations are Image reconstruction 
and Smoothing/Filtering. 

+ Image reconstruction. 

For a set of image points, usually grid points, deter- 
mine the image value, formally the average field value 
within the corresponding gridcell. In practice a few dif- 
ferent strategies may be followed 

- Formal geometric approach 

- Monte Carlo approach 

- Singular interpolation approach 

The choice of strategy is mainly dictated by accuracy 
requirements. For WVF we use the Monte Carlo approach 
in which the grid density value is the average of the 
DTFE field values at a number of randomly sampled 
points within the grid cell. 

+ Smoothing and Filtering: 

A range of filtering operations is conceivable. Two of rel- 
evance to WVF are: 

- Linear filtering of the field / 

Convolution of the field f with a filter function W a (x, y), 
usually user-specified, 



- Natural Neighbour Rank-Ordered filtering 
(see sec. I3.2|l . 

• Post-processing. 

The real potential of DTFE fields may be found in sophisti- 
cated applications, tuned towards uncovering characteristics 
of the reconstructed fields. An important aspect of this in- 
volves the analysis of structures in the density field. The 
WVF formalism developed in this study is an obvious ex- 
ample. 




(D5) 
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